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Abstract. - We performed molecular dynamics simulations to investigate the clustering 
instability of a freely cooling dilute gas of inelastically colliding disks in a quasi-one-dimensional 
setting. We observe that, as the gas cools, the shear stress becomes negligibly small, and the 
gas flows by inertia only. Finite-time singularities, intrinsic in such a flow, are arrested only 
when close-packed clusters are formed. We observe that the late-time dynamics of this system 
are describable by the Burgers equation with vanishing viscosity, and predict the long-time 
coarsening behavior. 



Introduction. - A gas of hard spheres is a standard model of statistical physics and kinetic 
theory [1]. It is surprising that a minor change in this model - the introduction of energy loss 
in the binary collisions - leads to consequences so dramatic. Among the many fascinating 
properties of the gas of inelastically colliding hard spheres [2], the clustering instability [3,4] 
plays a special role. No matter how small (but finite) is the inelasticity of collisions, the 
homogeneous cooling state (HCS) of this gas is always unstable if the system size L is large 
enough. When studying macroscopic properties of matter, a physicist deals, first of all, with 
the thermodynamic limit L — > oo. From this perspective the gas of elastically colliding hard 
spheres is a singular limit of the inelastic gas problem. The practical importance of the 
inelastic gas model stems from its being the simplest model of granular flow [2,5,6]. 

The clustering instability of a freely cooling inelastic gas involves formation of clusters of 
particles and generation of vortices [3,4,7]. The basic physics of the initial stage of cluster 
formation is simple: the inelastic cooling of the gas causes a pressure drop in the regions of 
enhanced density. This pressure drop drives an inflow of gas from the periphery and therefore 
provides a positive feedback to the instability. A traditional framework for quantitative theory 
here is granular hydrodynamics (GH), which assumes scale separation and is derivable from 
the Boltzmann equation, properly modified to account for the inelasticity of collisions [2]. 
Though the general criteria of its validity remain controversial [8], GH is well established 
at least when the following two criteria are met: (i) the granular gas is dilute, na D <C 1, 
and (ii) the particle collisions are nearly elastic, q <C 1. Here n is the local number density 
of the gas, a is the particle diameter, D > 1 is the dimension of space, q = (1 — r)/2 is 
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the inelasticity of collisions, and r is the coefficient of normal restitution (assumed constant 
throughout this paper). In this case the GH equations, linearized around the HCS, provide 
an accurate theory of the initial stage of the instability in terms of the (linear) hydrodynamic 
modes of the system: the shear mode (spontaneous formation of vortices) and the entropy, or 
clustering mode (formation of clusters) [3,4,7]. 

The further evolution of the instability is a hard problem. One difficulty here is technical, as 
the growing nonlinear shear and clustering modes become strongly coupled. Another difficulty 
is conceptual: GH breaks down in high-density regions. All previous attempts to develop a 
theory beyond linearization of the GH equations attempted to circumvent these difficulties. 
Bcn-Naim et al. [9] considered point-like particles inelastically colliding on a line. This strictly 
one-dimensional (ID) geometry makes a GH description problematic [8]. Still, Ben-Nairn et 
al. observed that, at long times, the ID system is describable by the Burgers equation with 
vanishing viscosity. Ernst et al. [10] considered a small two-dimensional (2D) system, where 
the entropy mode is suppressed, and dealt with the unstable shear mode. Baldassarri et 
al. [11] also studied instability of the velocity field in a homogeneous gas, by introducing a 
lattice model. 

Focusing on the entropy, or clustering, mode, Efrati et al. [12] had the shear mode sup- 
pressed, by working with a quasi-lD setting, see below. They solved numerically the low- 
density GH equations and found that, as the unstable system cools down, the shear stress 
becomes negligibly small, and a flow by inertia sets in. Formally, this flow develops a finite- 
time singularity: the velocity gradient and the gas density diverge at some location. Efrati et 
al. argued that the flow by inertia is an important intermediate stage of the clustering insta- 
bility. However, they did not address still later stages of the instability, when finite-density 
effects come into play. We report here the first MD simulations which probe the quasi-lD 
clustering of a dilute inelastic gas. We find that the simulated low-density stage of the clus- 
tering instability is in excellent agreement with the hydrodynamic predictions of Ref. [12]. 
What happens at later times? We observe that the attempted singularities are arrested only 
when hexagonally packed clusters are formed. The still later dynamics are describable by 
the Burgers equation with vanishing viscosity, in a striking analogy with the purely ID result 
by Bcn-Naim et al. [9]. Based on this observation, we predict the very-late-time coarsening 
dynamics of the system, all the way to its simple final steady state. 

Model system and clustering instability. - Consider an assembly of N identical hard disks 
of mass to = 1, diameter u = 1, and inelasticity < q <C 1 in a 2D box with dimensions 
L x and L y (L x 3> L y ). The initial number density of the gas no = N/(L x L y ) is very small 
compared to the hexagonal close packing density n c — 2/(v / 3c 2 ) — 1.155. The initial particle 
velocity distribution is gaussian with temperature To = 1. The parameters are chosen so that 
the low-density GH equations [5] are accurate until relatively late times, when the local value 
of no 2 is not a small parameter anymore. The boundary conditions are periodic at x = 
and x = L x and elastic (specular reflection) at y = and y = L y . Assuming a homogeneous 
cooling, one arrives at Haff's law [5] T(t) = T (l + t/t )- 2 , where t = (27r 1 / 2 crn gT 1/2 )- 1 
is the characteristic cooling time. The GH equations, linearized around the HCS, show that 
the HCS is unstable if L x and/or L y are large enough [3,4,7]. The density perturbations 
grow in time algebraically. The temperature and velocity perturbations decay, but the decay 
rates are smaller than that described by Haff's law. As a result, the flow tends to become 
supersonic [7]. A strong quasi-lD instability requires two criteria: an^Lxq 1 !" 1 3> 27T 1 / 2 , and 
anoLyq 1 / 2 <C 7T 1 / 2 ; the latter one guarantees suppression of the shear and clustering modes 
in the y-direction [3,4,7, 12]. 
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Fig. 1 - The local area fraction v(x,t) = n(x,t)/n c at times (a), 274 622 (b), 549 554 (c), 771055 
(d), 2.14972 x 10 6 (e), and 2.53252 x 10 6 (f). 



Flow by inertia. - We performed extensive event-driven MD simulations in this regime of 
parameters. We verified that no structure in the y-direction appears, as expected. Therefore, 
our diagnostics focused on ID coarse-grained fields: the density n(x,t), the mean velocity 
v(x,t) = (v x ,v y ) and the x- and y- components of the velocity fluctuations: T x (x,t) and 
T y (x,t). We report here a typical simulation with N = 12 500, q = 0.04, L x — 5 x 10 5 and 
L y = 25. The initial gas density is no = 10~ 3 , while to ~ 7.05 x 10 4 . Strong clustering 
instability is clearly seen in Fig. Q which exhibits formation of multiple clusters, and Fig. 
12 which shows an inflow of gas into the forming clusters. Meanwhile, the gas temperature 
T = T x + T y (not shown) rapidly decays with time. Throughout the simulated dynamics T x 
and T y remain close to each other, while their spatial inhomogeneity is small compared to the 
strong inhomogeneity of v x and n. 

The average particle energy versus time, E to t(t), follows Haff's law at early times, but 
deviates from it at later times, when strong hydrodynamic motions develop, see Fig. [21 (left) . 
In the hydrodynamic description E to t(t) — (tiqLx)^ 1 J q Lx (Et + E mac )dx, where Et = nT 
is the thermal energy density, and E mac — nv 2 /2 = n(v 2 + v 2 )/2 is the macroscopic kinetic 
energy density. The role of each term is elucidated in Fig. |31 Both Ex{t), and E mac (t) 
initially decay; Ex(t) decays faster. At later times Exit) and the y-component of E mac (t) 
continue to decay rapidly, while the x-componcnt of E rnac (t) decays much slower. As a result, 
Etot(t) is dominated by Et at early times and by the ^-component of E mac (t) at later times. 
Remarkably, Exit) continues to follow Haff's law quite closely until the latest simulated times. 

Figure 0Jl shows the time history of a typical cluster (the leftmost density peak in Fig. 
Q]c-f). The rapid density growth is saturated when n max approaches n c . A snapshot of 
the density peak region (Fig. [SJl indeed shows almost perfect hexagonal packing. The rapid 
density growth is shown in more detail in Fig. ^> which depicts 1 /n max versus time. The linear 
dependence, observed at intermediate times, indicates an "attempted" finite-time singularity 
nmax ~ (const— t)^ 1 . The same density singularity was observed in hydrodynamic simulations 
[12]; it is caused by a flow by inertia which develops when the forces acting on a fluid element 
vanish. The reason for it in the freely evolving inelastic gas is the continued rapid cooling, 



4 



EUROPHYSICS LETTERS 



v x 




0.05 - 



-0.05 - 



le+05 2e+05 3e+05 4e+05 5e+()5 le+05 2e+05 3e+05 4e+05 5e+05 




0.05 



-0.05 



/ " (c) - 



0.05 



-0.05 









(d) - 



le+05 2e+05 3e+05 4e+05 5e+05 le+05 2e+05 3e+05 4e+05 5e+05 



Fig. 2 - The horizontal velocity profiles at times (a), 274622 (b), 771055 (c), and 2.53252 x 10 6 
(d). The straight lines in Figs, c and d have a slope l/(t + C), where C = 3.8 x 10 . 



which makes the pressure and viscous stresses negligible [12]. The flow-by-inertia equations 
read 

dv x dv x Qtl d(nv x ) 

— +^— = 0, (a) and _ + ^_ = 0.(b) (1) 
These equations are soluble analytically in Lagrangian coordinates [13, 14]: 



v x (x,t) = v (£) , (a) n(x,t) = 
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Fig. 3 - Left: the energy balance of the system versus time. E to t is the average kinetic energy of 
the particles, Et is the average thermal energy, and E mac is the average energy of the macroscopic 
motions. The Haff 's law is indicated by H . Right: The horizontal velocity profiles in a region around 
the leftmost peak. Shown are v x versus x (a), and v x versus x — (t + C) v x (x,t)(b) at times 294518, 
410692, 549554, 636990, 658371, 680869 and 701507 with C = 3.8 x 10 5 . 
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Fig. 4 - Evolution of the leftmost peak in Fig. Figure a shows i/ mol versus time. Figure b shows 
n max versus time and a linear fit. The inset shows the ^-coordinate of the peak versus time. 



where C is an arbitrary constant, Uq(£) = efoo(£)/d£, while i>o(£) and ?i-o(£) are the v x and 
n, respectively, at the "initial" moment of time t = —C. The relation between Eulerian 
coordinate x and Lagrangian coordinate £ is x = £ + «o(£) (t + C). The finite-time singularities 
of both the velocity gradient, and the density occur when the denominator in Eq. b 
vanishes for the first time; it requires «q(0 < 0. 

The results of our MD simulations quantitatively agree with the flow-by-inertia scenario. 
Using in Eq. (jSJa) the straight-line fit, shown in Fig. QJd, we obtain C ~ 3.8 x 10 5 and 
|fo(£*)l — 1.0 x 10~ 6 , where is Lagrangian coordinate of the density peak. Now we verify 
Eq. © by observing [see Fig. |21 (right)] collapse of v x (x, t) versus £ = x — (t + C) v x (x, t), at 
different times, in the region of the leftmost density peak. The slope of the obtained straight 
line yields an independent estimate of |«o(£*)| which agrees within 3 percent with the value 
of 1.0 x 10~ 6 found earlier. Still another prediction deals with the shape of the density peak 
close to the time of attempted singularity. The flow-by-inertia theory predicts a density profile 



-2/3 



[14]. The simulations confirm this prediction, see Fig. [S] (right). 



As we observed, the density growth in the clusters is suppressed only when the density 
approaches n c . At a fixed temperature, the pressure of an assembly of nearly elastic hard 
spheres diverges like (n c — n) _1 as n approaches n c . Is there a pressure "revival" at n — > n c ? 
We used the empiric relation p = nT(n c + n)/ (n c — n), which interpolates between the dilute 
limit and the close-packing limit [15], to calculate the pressure field p(x,t). We found that, 
though Umax approaches n c , p{x,t) continues to decay with time, apparently because of the 
very rapid temperature decay. Therefore, the pressure revival mechanism can be ruled out. 

The Burgers model. - As the system approaches close packing it gets jammed. Faster 
particles cannot overrun slower ones, and the singular growth of the velocity gradient and 
density is arrested. A natural continuum model for a jammed flow is the Burgers equation [13] 
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Fig. 5 - Left: a snapshot of the system in the region of the leftmost peak at time 771 055. Right: 
The right tail of the density peak close to the time of attempted singularity. Shown is the log-log plot 
of the density profile n of the leftmost density peak versus x — x max at time 732 951. The dashed line 
is a power law 0.2 x [x — i m „)~ 2 ' 3 . 



in the limit of vanishing viscosity v — > 0, together with the continuity equation £[Jb. We 
shall call this model "the Burgers model" . Equation © is soluble exactly by the Hopf-Cole 
transformation [16,17]. The zero- viscosity limit of the solution, which is the subject of our 
interest, has the following properties. The solution is identical to that predicted by the flow- 
by-inertia model until the time moment when the flow by inertia would have developed a 
singularity. The attempted singularity gives way, in the Burgers model, to a "shock" : a jump 
in the velocity field which carries a density peak (cluster). In this coarse-grained description, 
the close-packed clusters have zero sizes but finite masses. At long times, when the shocks 
have "matured" [17], the quantitative predictions of the Burgers model are especially simple 
and can be conveniently tested in our MD simulations. One prediction is that dv x /dx is equal 
to l/(t + C) everywhere between the shocks, where C is the same constant (~ 3.8 x 10 5 in the 
reported simulation) as above. This prediction is tested in Fig. [2t and d. While the agreement 
in Fig. c is only fair (as the shocks have not yet matured), it improves considerably in Fig. 
d. Another prediction that we verified (see the inset of Fig. 0Ji) is that each shock moves 
with a constant speed until it collides with another shock. Furthermore, a collision between 
two shocks leads to their merger, without any thermalization of the system. Such a merger 
event can be seen in the left part of Figs. H c an d f- Finally, we verified that the zeros of v x , 
belonging to the "ramps", stay at rest as expected [17]. 

Late-time coarsening dynamics. - Event-driven MD simulations become very slow once 
the hexagonal close packing in the clusters is achieved. Fortunately, the times reached in 
our simulations were large enough to have verified the Burgers model as a proper late-time 
continuum model of the quasi-lD clustering process. Therefore, we can give a detailed pre- 
diction of the still later coarsening dynamics of the system, without a need to simulate it. If 
the number of clusters is large, the coarsening dynamics, which proceed via cluster mergers, 
can be addressed statistically. At this level of description the problem coincides with that 
of the decaying Burgers turbulence [16,17], or the ballistic agglomeration model [18]. For 
uncorrelated initial conditions, the average cluster mass grows with time like t 2 ^ 3 , the average 
velocity decreases like i" 1 / 3 , and the average distance between two neighboring clusters grows 
like i 2 / 3 [18]. This yields a long-time asymptotic scaling law E tot ~ E Vx ~ i~ 2 / 3 . At this 
stage, E tot is dominated by the kinetic energy of the clusters. The energy decay stops when 
an ultimate steady state of the system is reached: a single "super-cluster", moving with a 
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small constant speed (which can be determined from the initial data, by employing the mo- 
mentum conservation in the x-direction and assuming that all N particles are absorbed by 
this super-cluster). 

Summary and Discussion. - Our MD simulations fully support the hydrodynamic flow- 
by-inertia scenario [12]. The attempted singularities of this flow are suppressed when almost 
perfect hexagonally packed clusters are formed. At still later times the dynamics are describ- 
able by the Burgers equation with vanishing viscosity, in a striking analogy with the purely 
ID results by Ben-Nairn et al. [9]. Therefore, clustering instability in a quasi- ID setting is now 
well understood. What about a fully multi-dimensional geometry? It has been conjectured 
that, prior to the first attempted singularity the system should be describable by a multi- 
dimensional flow-by-inertia model [12]. One can expect that an interplay between vorticity, 
produced during the early stage of the instability and jamming provides a saturation mech- 
anism for the multi-dimensional singularities. Unfortunately, the multi-dimensional Burgers 
model [16, 17] assumes a potential velocity field and misses the important physics, caused by 
the presence of vorticity. The formulation of a continuum model free of this flaw should be 
the next step of theory. 

* * * 
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